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We investigate the exact solution of the honeycomb model proposed by Kitaev and derive an 
explicit formula for the projector onto the physical subspace. The physical states are simply char- 
acterized by the parity of the total occupation of the fermionic eigenmodes. We consider a general 
lattice on a torus and show that the physical fermion parity depends in a nontrivial way on the 
vortex configuration and the choice of boundary conditions. In the vortex- free case with a constant 
gauge field we are able to obtain an analytical expression of the parity. For a general configuration 
of the gauge field the parity can be easily evaluated numerically, which allows the exact diagonal- 
ization of large spin models. We consider physically relevant quantities, as in particular the vortex 
energies, and show that their true value and associated states can be substantially different from 
the one calculated in the unprojected space, even in the thermodynamic limit. 
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I. INTRODUCTION 

The Kitaev honeycomb model, with several variations, 
has attracted a lot of attention over the last years 
Many different interesting aspects of it have been studied 
in detail in the original work of Kitaev.^ There, an exact 
method of solution of the model based on the mapping to 
Majorana fermion operators is discussed. Although alter- 
native mappings and approximation techniques exist, ^ ^ 
Kitaev's method is widely applied, being ideally suited 
to this class of spin models. Further, the presence of 
an abelian and a non-abelian phase (in the presence of 
an external magnetic perturbation) was demonstrated.^ 
The Kitaev honeycomb model has a wide spectrum of 
physical applications, ranging from the description of 
strongly correlated materials^ to the analytical study of 
critical quantum spin liquids.^ It is also of central im- 
portance in the context of quantum information theory 
since its gapped phase provides a perturbative realization 
of the toric code.*^ Extensions of the honeycomb model 
have been lately proposed as promising candidates for 
the realization of a topological quantum memory.!^ Al- 
though very challenging, its physical realization has b e- 
come closer to reality thanks to recent proposals .'^^'II' 

In this paper we examine the projection to the physi- 
cal subspace of the exact mapping to Majorana fermions 
proposed by Kitaev.^ As briefly discussed in Refs. 12-14, 
unprojected and projected models have different physical 
properties, especially the parity of fermions. We derive 
here for the first time an explicit and immediately ap- 
plicable representation of the projector in terms of the 
parity of physical fermions. As it turns out, the physical 
fermion parity depends in a nontrivial way on the config- 
uration of vortices and on the lattice topology. Applying 
the projection to specific cases, we find large differences 
between projected and unprojected physical quantities 
(e.g. ground state and vortex energies, or spin-spin cor- 
relation functions). Such discrepancies exist both in the 
gapped and gapless phase and can also survive the ther- 
modynamic limit. Our analysis is consequently essential 



for the exact numerical study via Kitaev's exact map- 
ping of large spin systems, especially if one wants to go 
beyond the small system sizes of about 20 — 100 spins cur- 
rently accessible to various numerical approaches such as 
direct diagonalizatiorP^H^ or density matrix renormal- 
ization group (DMRG).^ These numerical approaches 
become necessary if exact fermionization techniques are 
not applicable. 

The paper is organized as follows. In Sec. [n|we briefly 
review the honeycomb model and the exact mapping to 
Majorana fermions introduced by Kitaev.l^ In Sec. Ill 
we compute the parity of physical fermions with periodic 
boundary conditions and a generic vortex configuration, 
which represents our main result. Section |IV| contains 
some applications to specific cases and Sec. |V|our final 
remarks. 



II. MODEL AND EXACT MAPPING 

The Kitaev honeycomb model is a quantum compass 
model:^ defined on an hexagonal lattice A as follows 



H 



(1) 



where cr^ are the Pauli spin operators at site i e A 
{i = 1, . . . , 2A^). In Eq. Q, the sum runs over all the 
pairs of nearest-neighbor sites and the directions of the 
Ising interactions are determined by the orientations of 
the corresponding links {aij = x^y^z for x-, y-, z-links 
respectively, see Fig. [T]) . 

To solve this spin model in an extended Hilbert space 
£, one can associate at each site i four Majorana modes 
Ci^bf ^b^ } By defining = ibfci^ the original Hamil- 
tonian in Eq. ([T]) is mapped to 



(2) 
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FIG. 1. Honeycomb lattice, with basis vectors ei,2- The di- 
rections of X, z-Ymks> are indicated, as well as the vector 
T joining B and A sublattices (white and black dots, respec- 
tively) . The most general torus on the lattice can be specified 
by Liei and L2G2 + Mei (here Li = 4, L2 = M = 2). The 
numbers 1,2, ... 16 label the sites as described in the main 
text, see Eq. ([s}. 

where for nearest-neighbor sites Aij = JaijUij and 

These operators satisfy u^j = 1. Furthermore, they all 
commute with each other and also with H. Therefore, 
the extended Hilbert space splits into C = ®uCui where 
u represents a configuration of Uij = ±1. Notice that 
Uij = —Uji. So, whenever we specify the values of Uij^ we 
assume conventionally that i is in the A sublattice (see 
Fig. [1]). In each subspace Cu^ the operator matrix Aij 
are replaced by numbers A^j and Eq. ^ thus describes 
non-interacting Majorana fermions. 

The eigenmodes can be easily obtained with a canoni- 
cal transformation to new Majorana operators 

(61, 61,. ..,6^, 6^) = (ci,...,C2Ar)Q'', (4) 

which for the specific configuration u brings H to the 
form Hu = ^^rn^rn^mPm-i whcrc 6^ are the positive 
eigenvalues of 2^74^. By introducing the fermion opera- 
tors am = 1/2(6^ + ib^) and rim = cim^m we obtain 

H^ = ^em{nm- 1/2), (5) 

m 

with ground state energy £^0 = "l/^Xl^e^. The or- 
thogonal matrix will have a crucial role in the follow- 
ing to obtain the projection operator. 

III. PHYSICAL FERMION PARITY 

The key advantage of Kitaev's solution is to reduce the 
problem of finding the eigenvalues of a 2^^ x 2^^ ma- 
trix to the diagonalization of the 2N x 2N matrices A'^. 



However, the final spectrum and eigenstates are in the ex- 
tended Hilbert space £, and a projection V to the phys- 
ical subspace is necessary.^D The physical states satisfy 
Dil"^) = 1^) for all the gauge operators Di = bfb^b^Ci 
and the explicit form of V is^ 

^=n(H^)=iEnA, (6) 

where the summation runs over all possible subsets of 
indices {i\ C A. Within the physical subspace the a^'^'^ 
operators satisfy the usual algebra of Pauli matrices and 
therefore R and R are equivalent. 

To establish a more explicit formula for V one can 
note that, in the summation appearing in Eq. ([6|, the 
two terms corresponding to a subset {%] and its comple- 
mentary set A \ {%) simply diff er by a factor Hii^i ^i- 
Therefore, V factorizes as follow4^^^ 

(7) 

where the prime indicates that the summation in S (in 
the first parenthesis) is restricted to half of all possible 
subset of indices: if {%] is included, then A \ {%) is not. 

We then consider Hi^^i projector T^o [the sec- 

ond parenthesis of Eq. ([7|)]. From the definition of Di, it 
clearly consists of a product of all the q and 6^'^'^ op- 
erators. By applying the anticommutation rules we can 
pair corresponding 6^'^'^ operators, and express them in 
terms of the conserved quantities Uij. To do this, it is 
necessary to know the topology of the lattice, from which 
the correct pairing is determined. We consider here a 
model defined on a torus with basis vectors Liei and 
^262 + Mei, as illustrated in Fig. [l] for a special case 
(Li = 4, L2 = M = 2). This represents the most general 
choice of periodic boundary conditions and TV = L1L2. 
It is also necessary to fix the correspondence between 
i = 1, ...,2N and lattice sites. We fix the labeling as in 
Fig. [1] By taking the origin on site i = 1 (on the B sub- 
lattice) the position of the sites with odd values of i is 
given by 

where {a\b) indicates the integer division and (a mod b) 
the reminder. For even the position is = r^_i + r 
(see Fig. [l]). After pairing the 6^'^'^ operators into the 
Uij , the result is proportional to the parity operator tTc = 
{—i)^ Yli^i' We then express this quantity in terms of 
the eigenmodes: 

TTc = det(Q^)7r, (9) 

where tt = nm=i(-'- ~ 2nm) is the parity of the eigen- 
modes am- A proof of Eq. (|9| is provided in Appendix [A| 
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Finally, we find for Vq, 
2Po = 1 + (- 



-l)^det(Q«)7rJJ' 



(10) 



where 6> = Li + L2 + M{Li — M) and tt has eigenval- 
ues +1(— 1) if the total number of physical fermions is 
even (odd). A complete derivation of Eq. (10) is given 
in Appendix [BJ It is important to notice that, in apply- 
ing Eq. (10), the labeling of the lattice described above 



should be used. For example, det(Q^) depends on the 
choice of the labeling. 

We would like to point out now the differences between 
Eq. (10) and other discussions in the literature.l^^lt^ 



Firstly, the parity of physical fermions (the only relevant 
ones) TT is the parity of the eigenmodes am and not simply 
TTc which is of no calculational use. Eq. (|9| shows now the 
precise relation between tt and tTcI for a certain configu- 
ration the two parities are different if det((5^) = — 1. 
We note also that tTc is not a gauge invariant quantity, 
while the physical parity tt obviously is, i.e., [ii^Di] = 0. 
Therefore Eq. ([9| allows one to form the gauge invariant 
quantity det((5^ f^^- ^-^ ii^j. A second feature revealed 
by our analysis is that tt (and tTc as well) depends in a 
nontrivial way on the boundary conditions through the 
factor (—1)^, which does not appear in Refs. 14 and 15, 
Being directly applicable to the eigenstates G >Cn, 
identified by their occupation numbers (1 — 2n^) = ±1, 
Eq. ([To| is extremely convenient: it immediately shows 
whether Vo gives or 1 on the eigenstate |^)n- In the 
former case, the state is clearly unphysical. In the sec- 
ond case, V\'^)u = S\'^)u 7^ since, as seen in Eq. ([7|), 
the 2^^~^ terms of S all correspond to different con- 
figurations of Uij. In conclusion, Eq. (10) is sufficient 
to determine if 1^)^^ has some overlap with the physical 
subspace or lies completely outside of it, and makes clear 
that the crucial quantity is the parity of physical fermions 
tt: physical states have either even or odd occupation of 
the eigenmodes am depending on both the configuration 
u and the choice of boundary condition. 



IV. EXAMPLES OF PROJECTED STATES AND 
ENERGIES 

We discuss in this section a few examples illustrating 
the difference between physical and unphysical results. 
These examples should make clear that it is not sufficient 
to calculate the properties of the system in the extended 
space, but the projection must be carefully applied to 
take advantage of all the power of Kitaev's exact map- 
ping. While we focus here on the ground state and vortex 
excitation energy, we expect that similar discrepancies 
exist for other physical quantities. 

By applying Eq. (10), the only factor which is not im- 



mediately found is det(Q^) and we show in the following 
how it can be explicitly evaluated when Uij = 1. The fi- 



nal result, Eq. (18), nicely complements Eq. (10) for this 



vortex- free sector. 




FIG. 2. An illustration of the partitioning into Q, Q± de- 
scribed in the main text. The vectors bi,2 define the recipro- 
cal lattice and the hexagon is the first Brillouin zone. The four 
square points are wavevectors in Q while the dots and crosses 
are in and Q- , respectively. In this example Li = L2 = 6 
and M = 0. 



For an arbitrary configuration of the Uij^ the Fourier 
transformation cannot be used to calculate analytic re- 
sults. However, det((5^) can be determined numerically 
with negligible computational effort. This allows us to 
obtain the exact numerical solution of the spin Hamil- 
tonian at very large values TV and to explore the effect 
of the projection when approaching the thermodynamic 
limit. 



A. Vortex- free sector 

In the vortex-free sector with Uij = 1 we can proceed 
as in Ref. [T] by making use of the Fourier transform on 
the lattice 




(11) 



where the positions of the B lattice sites are given in 
Eq. ([8|. We now consider a partition of the N possible 
values of q (fixed by the periodicity of the lattice) in three 
sets: Q and Q±. We say that q G ^1 if ±q are the same 
(up to reciprocal lattice vectors), ft contains at most four 
wave vectors, depending on L 1^2 and M, and always con- 
tains 0. The remaining wavevectors can be partitioned 
in a way that ±q always belong to two distinct sets Q±. 
An example of such partitioning is illustrated in Fig. [2] 
for the special case Li = L2 =6 and M = 0. We define 
new Major ana modes as 7^ = A/2aqA for q G and 



7^1 



-qA? 



7q,2 =^(^^qA-0^-qA), (12) 



for q G where X = B refers to the sublattice 
(see Fig. [T]). This canonical transformation of the q 
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can be constructed in two steps. First we rearrange the 
(ci,C2, . . . ,C2Ar) into 



(C2, C4, . . . , C2N-, Cl, C3, . . . , C2Ar-l), 



(13) 



a transformation which has determinant ( — 1)^(^+1) Z^. 
Notice that, in Eq. (13), the modes are partitioned be- 



tween A (first half) and B (second half). Furthermore, 
for each sublattice the same order of z-dimers appears, 
since q and q+i (with odd i) belong to the same z-dimer 
(see Fig. [1]). Because of this structure, the second trans- 
formation [from Eq. (13) to the 7q,7q^Q, modes] has two 
identical blocks labeleaoy A = A, B and the determinant 
is simply 1. 

The Hamiltonian, rewritten in terms of the new Majo- 
rana modes, is diagonal in q and its coefficients are given 



by /(q) = 2(J, 



Jy 



Jz) and its complex 



conjugate. A further diagonalization with respect to the 
index a of is achieved with the rotation of the B 
operators: 



sm 



1) COs(0q)y I 7 



7.^ 



'q,2 



(14) 



where is the phase of /(q), i.e., /(q) = |/(q)|e*'^'*. 
This transformation has again determinant 1 and brings 
Hu to 




l/(q)l 



E 

a=l,2 



1^ 7^ 

/q,Q; /q,Q; 



(15) 

Finally, as discussed below Eq. (|4|, brings the Hamil- 
tonian to the form Hy = | Xlg €(q)^q^q with e(q) > 0. 
This can be achieved in Eq. (fTsj) by relabeling the Majo- 
rana operators. If q G 1^+: 



If q 



7q,l,^q 
7q,2,C 



7^' 
= 75- 



(16) 
(17) 



BiA) 
7q 



7q 7 



e n and /(q) > 0(< 0): .q - ;q , .q 

. The determinant of this last transformation is 
(^ — l^x-\-NiN-i)/2^ where x is the number of reciprocal lat- 
tice vectors q G 1^ such that /(q) < 0. By combining this 
factor with the one from Eq. (13) we obtain 
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FIG. 3. Physical (solid line) and unphysical (dashed line) fi- 
nite size corrections to the thermodynamic energy of a vortex- 
free configuration Uij — 1 for Li — L2 — L. The main plot 
a) refers to the gapless phase with M = 0, Jx — Jy — Jz — 
and £0 — —1.5746.^ The inset h) corresponds to the gapped 
phase with M = 1, J^; = = 0.2, J^ 1, and £0 - -1.0202. 



is the energy per unit cell in the thermodynamic limit. 
We consider in Fig. [3] a square lattice (Li = L2 = L) 
with Uij = 1 and two different choices of boundary con- 
ditions and couplings, and plot SE{N) as a function of L 
{N = in this case). The original result calculated in 
Ref. T is reproduced in the main panel of Fig. [3] (dashed 
lines) and evidently refers to the unphysical energy, which 
always underestimates the correct result. The physi- 
cal and unphysical energies are always distinct, unless 
= for some fermion mode. Being in the gapless 
phase, the difference approaches zero at large system size 
as 1/L. The inset represents an example in the gapped 
phase: remarkably, since there is always one fermion in 
the physical ground state, the difference between pro- 
jected (solid line) and unprojected (dashed line) results 
survive the thermodynamic limit where the true energy 
correction does not vanish. That the state with zero 
fermions is never physical for any L is immediate from 
our analytic result Eq. (18) since x = iii the gapped 



phase [/(q) > 0]. On the other hand, projected and un- 
projected states are different (have different parity) in 
the thermodynamic limit in both the gapless and gapped 
phases. Therefore our projection protocol is necessary 
to determine the physical quantities of the model in both 
the gapless and gapped phase even in the thermodynamic 
limit. 



det (Q") = (-1) 



for 



1. 



(18) 



Notice that x depends in a non trivial way on the bound- 
ary conditions Li^2,^, and the couplings Jx,y,z- Never- 
theless, for a given choice of the model, it can be easily 
computed. 

Following Ref. 1 we examine now the finite size correc- 
tion SE{N) = Eo{N) — sqN to the ground-state energy 
Eo{N) with Uij = 1, where 



So = lim Eo{N)/N 

N^oo 



(19) 



B. Energy of two adjacent vortices 

We consider next the energy to create vortices in the 
system. These are present on hexagonal cells for which 
the product of the six Uij is —1. As an interesting ex- 
ample we study configurations with two adjacent vor- 
tices, obtained by setting Uij = —1 for a single link. The 
ground state energy of such two-vortex configuration of 
the Uij is denoted as E2{N) while, as before, Eo{N) is the 
ground energy with all Uij = 1. We define the excitation 
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FIG. 4. Physical (thick soUd hne) and unphysical (thick 
dashed hne) excitation energy of two adjacent vortices. This 
plot refers to the gap less phase with J^, = = = 1, 
Li =2L, L2 — M — and a single Uij — —1 with ij being a 
z link. Thin dashed lines extrapolate to the thermodynamic 
limit. The dots at L = 3 (unphysical) and L = 13 (physical) 
represent the largest system sizes with negative excitation en- 
ergy. 



energy of a vortex as 



ilE{N) = -[E2{N) - Eo{N% 



(20) 



which is plotted in Figs. |4]and|5]as function of L for dif- 
ferent choices of the parameters (given in the captions). 
In particular, Fig.|4]refers to the gapless phase and Fig.|5] 
to the gapped phase. Notice also that it is important in 
general to specify if Uij = — 1 refers to an x, or z link. 

Since the vortex state is not translationally invariant, 
we obtain the energy spectrum and det((5n) numerically. 
This can be done efficiently, since it only involves 2Nx2N 
matrices. In Figs. |4] and [5] we show physical excitation 
energies up to L = 26 and L = 35, respectively (corre- 
sponding to = 2704 and 2L^ = 2450 physical spins), 
obtained with a standard tabletop computer and high 
level language routines (matlab). We have checked that 
for the lowest possible system size the spectrum obtained 
with this method is in agreement with direct diagonal- 
ization of the physical model. However, direct diago- 
nalization is limited to systems with only a few tens of 
spins, by making use of ver y intensive parallel comput- 
ing numerical routines.l^^^ Other numerical approaches 
such as DMRG^ allow to address larger spin systems 
(2A^ < 100) but still much smaller than the ones accessi- 
ble with the projection protocol presented here. 

Similarly to Fig. [3j physical and unphysical results are 
generally different. In Fig. |4] (in the gapless phase) they 
become equal only if there is a zero energy mode or in the 
thermodynamic limit. On the other hand, finite size cor- 
rections are very important: AE{N) shows pronounced 
oscillations with an amplitude which is of the same order 
of magnitude of A£^(oo), for up to a few thousand spins. 
Remarkably, such oscillations result in negative excita- 
tion energies of the vortex pair up to 676 spins (instead 
of 36, in the extended space). 

In the gapped phase, the difference between projected 
and unprojected results can survive the thermodynamic 
limit. In Fig. [5^ the oscillations in the physical solution 
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FIG. 5. Physical (solid line) and unphysical (dashed line) 
excitation energy of two adjacent vortices in the gapped phase 



with Jx = Jy = 0.1, Jz 

Uij = — 1 with ij being s 
while panel 6) to M = 1. 



= 1, Li = L2 = 
y link. Panel a] 



-- L, and a single 
refers to M = 



(solid line) persist for L ^ 00 and their size is equal to 
half the gap of the fermions (~ J^), much larger than 
the excitation energy of the vortex pair in the unphysical 
space (dashed line). A different choice of boundary con- 
ditions can lead to a situation where is always large 
and negative (~ —Jz) and has a well defined thermody- 
namic limit, as illustrated in Fig. ^jp. 

The thermodynamic limit of the unphysical energies 
(dashed curves) of Fig. [5] is well approximated by a high 
order perturbative expansion, which allows to derive an 
accurate effectiv e Hami ltonian including vortex energies 
and interactions l "*^ * ^ * ^ * ^ However such low-energy Hamil- 
tonian does not contain explicitly the constraints on the 
allowed vortex configurations, and some care is necessary 
to establish which states are physical. For example, the 
vortex-free configuration of Fig. [sJd does not belong to 
the low-energy subspace and the two- vortex state has al- 
ways lower energy. It is worth pointing out here that the 
negative vortex energies encountered in these examples 
are not in contradiction with Lieb's theorem as originally 
formulated in Ref. |27| (see also Ref. l28p. 



V. CONCLUSION 

We have obtained here an explicit form of the pro- 
jection operator which allows us to extract the physical 
properties of the honeycomb model for large lattices. The 
parity of fermions in the physical sector is directly given 
through our Eq. (10) and depends in a nontrivial way 
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on the vortex configuration and tiie periodic boundary 
conditions. By applying Eq. (10), we have examined the 
energies of vortex-free and two- vortex configurations and 
showed that significant deviation from the physical val- 
ues can exist if the projection operator is not taken into 
account. Such differences between projected and unpro- 
jected quantities can persist up to large values of N and 
they can survive the thermodynamic limit. 

Applying the projection only requires to determine the 
parity of physical fermions. Therefore, it does not in- 
troduce any additional complication related to the sym- 
metrization over all gauge transformations ([6|, a proce- 
dure which never needs to be implemented in practice. As 
known, the energies of projected (i.e. symmetrized) and 
unprojected (i.e. unsymmetrized) states are the same 
and the spin correlation functions S^f{t) = (erf (t)cr^(O)) 
(whose certain exact properties were discussed in Refs. |2] 
and ii2Q]) can be conveniently computed with unprojected 
eigenstates |^)^i.'^This is possible thanks to the fact that 
the spin operators are gauge-invariant. However, one 
should make sure that only states with 7^|'0)n 7^ are 
included in S^j^{t). Therefore deviations from the un- 
projected results exist for the S^j^ (t) as well. 

More generally it is obvious from our discussion that 
all dynamic and thermodynamic quantities derived from 
H (for example, the partition function), depend on the 
physical spectrum and eigenstates and thus differ from 
those of the unprojected model H. Therefore, we think 
that it would be interesting to apply our method toprob- 
lems which have been studied without project iorf^J^^^^ 
and compare the differences. Our work is generally rel- 
evant to spin models to which Kitaev's solution applies, 
like the honeyc omb model perturbed by a weak mag- 
netic interacting with cavity modes^ or with 
different link distributions^^, and a three-dimensional ex- 
tension of the honeycomb model recently proposed in 
Ref. |22] The case of open boundary conditions can 
also be simply obtained by extending Eq. ([T]) to site- 



dependent couplings Ja^j 
boundary) . 



Jij (and Jij = on the 
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sake of simpler notation we relabel {h'^^h'l^ ...^h'^^h'^) as 
(61, 62, ^2Ar). The relation between the c and h Majo- 
rana fermion operators is given by Eq. Q: 



Cj 



(Al) 



where Q is an orthogonal matrix. Let us consider the 
set S2N of the permutations of 1,...,2A/'. Since the q 
anticommute, we can write the product C1C2...C2N as a 
sum over all permutations a G S2N as follows: 



2N ^ 2N 



i=l 



(A2) 



k=l 



whe re e (a) is the sign of permutation a. By using 



Eq. (|A1| we can write 11^=1 <^cr(/e) as 



}2N,i2N ^icj{l)"'^icj{2N)-> 



(A3) 



where the numerical factor ]\^ Qk,ik only depends on the 
values of i/c, an d no t on the permutation a. This allows 
to express Eq. (A2) in the following form 



il,--,i2N k=l \ cr / 



(2iV) 



It is not difficult to check that the sum over a gives zero 
if the values of two of the indexes ik are equal. When the 
ik are all distinct, they are a permutation of 1,...,2A^: 
ik = o-\k) and i^r^k) = ^'(^(^))- Furthermore, we can 
anticommute the Major ana operators to the canonical 
order 6i62---^2Ar, which introduces the sign e((j)e(cr').This 
leads to 



M fe '^'^'^ n Qk,<r'ik)) hb2...b2N- (A5) 



(2iV) 



where the {2N)\ is canceled by the sum over a. The 
remaining factor is simply the determinant of Q, and the 
following relation is obtained: 



2N 



2N 



Y[ci = det{Q)l[b,, 



(A6) 



2=1 i=l 

from which Eq. ([9| directly follows. 



Appendix B: Derivation of Eq. (10) 



Appendix A: Derivation of Eq. ([9]) 

In this appendix we present a detailed derivation of 
Eq. (|9| which gives the relation between the parity tTc 
and the parity of the physical fermions tt. For the 



In this appendix we present details of the derivation of 
Eq. ([10|. From Eq. ^ we know that 



2N 



2Po = l + Y[Di 



(Bl) 



i=l 
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with gauge operators Di = bfb'^b^Ci. The product of 
ah gauge operators can then be rewritten in terms of 
Majorana operators as 



2N 



n A = bible . . . bl^bl^bl^C2N. (B2) 



where N = L1L2. In order to pair the b operators and 
form the corresponding u operators, we first move all the 
c operators to the right of the b operators. The fermionic 
anticommutation relations then lead to a phase factor of 
(-1)^1 with ^1 = ?>N{2N - 1): 



where = ^2- Note that the phase ^4 arises because of 
the boundary conditions along ei. To form the u\2Li op- 
erator (iifg in Fig.iri) one has to move 6f after , which 
introduces a (— 1) lactor. This procedure has to be re- 
peated L2 times, for all pairs of the form ^2nLi^2Li(n-i)+i 
with n = 1, 2, . . . , L2 (in Fig.jl] these are 6g6f and ^fg^g). 
Hence, the factor (—1)^^ arises. 

Finally, in order to identify 11^=1 ^ product of 

^ we first decompose it in L2 products of 2Li terms: 



2N 



(B7) 



2N 



Y[D, = {-lf^blb\bl 



2N 

blNblNblNWci. (B3) 



1=1 



Since 6f and b^j^^ are always separated by two fermionic 
operators, we can move all the b^ to the left without 
introducing any phase factor. We then group together 
all the b^ operators to the right of the 6^, at the cost 
of an additional phase (— l)'^^ where (j)2 = N{2N — 1) 
[notice that (-1)^1+^2 = 1]. 



2N 



2N 



2N 



2N 



2N 



Cl, 



(B4) 



k=l 1=1 



For the sake of clarity we explicitly write in this ap- 
pendix each u operator as where a refers to the link 
associated to u: Uij = u^-^ = ib^''^ b^'^ . In order to pair 
the 6" {a = x^y^z) operators into li", it is necessary 
to fix the correspondence between i = 1, ...,2A^ and the 
sites of the lattice. We choose here the labeling defined 
by Eq. (|8| and illustrated in Fig. [l] for a particular case 
(Li = 4 and L2 = M = 2). It is then straightforward to 
see that 



2N 
k=l 



(-1) 



03 



n 



(B5) 



is from the definition of ufj = ibfbj^ 



where the factor ' 

and (ps = N arises from the convention of specifying u^j 
with i e A. Therefore, the above have always the 
form u^' 



It is also not difficult to rearrange the b^ operators: 



where Tn 



n-=l&^L,(„-i)+, (e-g- 7I = 6?6l...6lfor 
Fig. fl]). As a first step, we rewrite each Tn by moving all 
the of with odd i on the left side (and keeping them in 
increasing order of i) while rearranging the b^ with even 
i (now on the right side of each Tn) in decreasing order. 
These operations do not introduce any additional phase 
factor in the final expressions. For example, in the case 
of Fig. [1] we can write Ti = blb^b^b^b^bQb^b2 . The ad- 
vantage of this ordering is that most pairs of b^ operators 
are now straightforward to form. For example, the pairs 
62^9, b\b\-^^ ... of Fig. [1] are now easily formed from the 
rearranged string ( |B7| ), of the form . . . bQb^b2bgb\-^b\^ .... 
The only difficulty is a remaining product of b^ 



-l)i^2N^2N- 



■2' ' ' ^2N- 



2(Li 



0, (B8) 



which requires some care in pairing to account of the toric 
boundary conditions. The rearrangement of Eq. (B8) 
introduces the phase factor (— 1)"^^ in the final expression: 



2N 



AN 



n 

{m,n) 



(B9) 



where ^5 M{Li - M) + Li. 

Collecting all the terms, we finally obtain 

^ (_l)03+04+^5 ^ 2Ar 



(BIO) 



i=l 



With the aid of Eq. ^ we then find 
2^0 = l + (-l)^i+^2+M(Li-M) det(Q^)^ Yl u,j, (Bll) 



2N 



(-1)^ 



AN 



n 



(B6) which is Eq. (10) of the main text. 
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